Theory and computation of covariant Lyapunov vectors 
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Lyapunov exponents are well-known characteristic numbers that describe growth rates of pertur- 
bations applied to a trajectory of a dynamical system in different state space directions. Covariant 
(or characteristic) Lyapunov vectors indicate these directions. Though the concept of these vectors 
has been known for a long time, they became practically computable only recently due to algorithms 
suggested by Ginelli et al. [Phys. Rev. Lett. 99, 2007, 130601] and by Wolfe and Samelson [Tellus 
59A, 2007, 355]. In view of the great interest in covariant Lyapunov vectors and their wide range of 
potential applications, in this article we summarize the available information related to Lyapunov 
vectors and provide a detailed explanation of both the theoretical basics and numerical algorithms. 
We introduce the notion of adjoint covariant Lyapunov vectors. The angles between these vectors 
and the original covariant vectors are norm-independent and can be considered as characteristic 
numbers. Moreover, we present and study in detail an improved approach for computing covariant 
Lyapunov vectors. Also we describe, how one can test for hyperbolicity of chaotic dynamics without 
explicitly computing covariant vectors. 

Keywords: covariant Lyapunov vectors; characteristic Lyapunov vectors; forward and backward Lyapunov 
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INTRODUCTION 

High-dimensional nonlinear systems like coupled oscil- 
lators, dynamical networks, or extended excitable me- 
dia often exhibit very complex dynamics that is diffi- 
cult to analyze and to characterize. From a practical 
point of view there only a few concepts have been de- 
veloped for studying low-dimensional systems that can 
efficiently be applied to high-dimensional attractors, too. 
An important example are Lyapunov exponents that de- 
scribe growth rates of perturbations applied to a trajec- 
tory in different state space directions. These exponents 
are a central point in the investigation of chaotic dynam- 
ical systems. They are related to a number of different 
physical properties such as sensitivity to initial condi- 
tions or local entropy production and can be used to es- 
timate the (Kaplan- Yorke) dimension of (even very high- 
dimensional) attractors 

Mathematically, Lyapunov exponents are defined in 
tangent space. This space is spanned by all possible in- 
finitesimal perturbations that can be applied to a state of 
the system. The dimension of the tangent space is equal 
to the dimension of the original phase space. In general, 
the tangent space is an inner product space, but often the 
tangent space is defined as an Euclidean space where the 
inner product is just the ordinary scalar (dot) product. 
The dynamics in this space is generated by linear opera- 
tors, that determine the evolution of perturbation vectors 
from one point on the trajectory to another. These oper- 
ators are called tangent linear propagators or resolvents. 
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The tangent space is a very important subject of study. 
On the one hand, the tangent space dynamics is closely 
related to the dynamics of the original system. One can 
obtain key characteristics of the original system observ- 
ing the associated tangent space dynamics. On the other 
hand, the tangent space is linear and the dynamics in 
this space is determined by the action of linear opera- 
tors. This means that analysis methods as well as results 
are universal for a wide class of systems. 

Besides the growth rates of perturbations the direc- 
tions of this growth are also important. There are var- 
ious concepts identifying these directions including bred 
vectors 0, , which are finite-amplitude perturbations 
initialized and periodically rescaled within the original 
phase space, singular or optimal vectors 0, which 
are the singular vectors of a finite-time propagator, or 
finite-time normal modes @, defined as eigenvectors of 
the propagator. 

Orthogonal sets of singular vectors related to the prop- 
agators operating on infinite time intervals were referred 
to by Legras and Vautard as forward and backward Lya- 
punov vectors Q • These vectors can be computed in par- 
allel with the Lyapunov exponents 0, @Ij and, thus, are 
closely related to them. Unlike the exponents, the for- 
ward and backward Lyapunov vectors depend on time, 
i.e., they are different for different trajectory points. An- 
alyzing the orientation of these vectors, one can expect 
to recover the local structure of an attractor. But un- 
fortunately, the forward and backward Lyapunov vectors 
provide only limited information. They always remain 
orthogonal and thus cannot indicate directions of sta- 
ble and unstable manifolds as well as their tangencies. 
These vectors are not invariant under time reversal and 
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are not covariant with the dynamics. The latter means 
that forward (or backward) vectors at a given point are 
not mapped by tangent propagators to the forward (back- 
ward) vectors at the image point. Another drawback of 
these vectors is their norm-dependence, i.e., they depend 
on the definition of the inner products and norms in the 
tangent space Q- 

The concept of norm- independent Lya.punov vectors ha 
been known for a long time [l|, 0, H, [13 ■ However, only 
recently two efficient algorithms for computing these vec- 
tors were suggested almost simultaneously by Wolfe and 
Samelson and by Ginelli et al. After Ginelli et al. 
we call these vectors covariant Lyapunov vectors. Note 
that these vectors are also referred to as characteristic 
Lyapunov vectors 0, [HI • These vectors are not orthogo- 
nal, they are invariant under time reversal and covariant 
with the dynamics in the sense that they may, in princi- 
ple, be computed once and then determined for all times 
using the tangent propagator. (Note that this is the case 
only for exact covariant vectors, while those computed 
numerically do not demonstrate perfect covariance due 
to the accumulation of numerical errors.) The covariant 
Lyapunov vectors can be considered as a generalization of 
the notion of "normal modes." They are reduced to Flo- 
quet vectors if the flow is time periodic and to stationary 
normal modes if the flow is stationary . 

In view of potential wide applications to the analysis of 
complex, high-dimensional dynamics, the covariant Lya- 
punov vectors receive a lot of interest of researchers |l3l ~ 
[l9fl . For these extensive studies to be productive, it is 
important to analyze the Lyapunov vectors systemati- 
cally. In this paper we summarize features of forward, 
backward and covariant Lyapunov vectors and provide a 
detailed explanation of both the theoretical basics and 
numerical algorithms. We present and study in detail an 
efficient method for computing covariant Lyapunov vec- 
tors, which can be considered as a modification of the 
method by Wolfe and Samelson. Moreover, our general 
approach reveals the existence of adjoint covariant Lya- 
punov vectors. This is not an independent type of char- 
acteristic vectors, because given the covariant vectors, 
one can always compute the adjoint ones. However, the 
angles between corresponding covariant and adjoint co- 
variant vectors provide a compact representation of the 
information contained in the covariant vectors and can 
be used as characteristic numbers. In particular, the 
presence of homoclinic tangencies is indicated by orthog- 
onality of corresponding original and adjoint covariant 
vectors. Since the covariant as well as the adjoint covari- 
ant vectors are norm-independent their angles also are 
invariant with respect to the norm. 

The structure of the article is as follows. In Sec. U 
we present the theory of Lyapunov exponents and for- 
ward and backward Lyapunov vectors, and in Sec. HIl we 
describe numerical methods for computing them. Sec- 
tion |llT] presents the theoretical aspects of covariant Lya- 
punov vectors, and in Sec. II VI we describe different meth- 
ods of computing covariant vectors. Finally, in Sec. |V] a 



simple illustrative example is presented. In Sec. IVII we 
summarize the results presented. 

I. LYAPUNOV EXPONENTS, FORWARD AND 
BACKWARD LYAPUNOV VECTORS 

A. Basic definitions 

Consider a system whose dynamics can be described 
by an ordinary differential equation 

u = g{u,t), (1) 

where u = u{t) Cz M'" is an m-dimensional state vec- 
tor that changes in time t, and g{u,t) G K"* is a non- 
linear vector function. We are primarily interested in 
high-dimensional systems, so m is assumed to be large. 
Equation ([1]) can model a system with many interact- 
ing point-wise dynamical elements, or it can be a finite 
step size approximation of a spatially extended system 
that appears after discretization of spatial derivatives. 
Infinitesimal perturbations to a trajectory of this system 
are described by the following equation: 

v^J{u,t)v, (2) 

where J{u,t) e ]^™x™ jg the Jacobian matrix composed 
of derivatives of the vector function g{u,t) with respect 
to components of the vector u. The fundamental matrix 
M e R"^^™ for Eq. (m can be found as a solution of the 
matrix equation 

M = J(M,t)M, (3) 

where any non-singular matrix can be used as an initial 
condition. 

The tangent linear propagator or resolvent is defined 

as 

j^{tiM) = mt2)mtir\ (4) 

and can be represented by a non-singular m x m matrix. 
The propagator evolves solutions of Eq. ^ from time ti 
to time t2'. 

v{t2)^ T{tiM)vit{), (5) 

where v{ti) and v{t2) are tangent vectors at times ti and 
^2, respectively, computed along the same trajectory of 
the base system ([1]). According to Eq. (j4|). the propa- 
gator is always non-singular and J^(ti,t2) = ^{t2,ti)~^ . 
Furthermore we define the adjoint tangent propagator: 

Q{tiM) = ^{tiM)'^ , (6) 

where "— T" denotes matrix inversion and transposition. 
In general, a non-Euclidean norm can be defined in the 
tangent space, so that instead of the transposition a gen- 
eralized adjoint with respect to the chosen norm has to 
be used. In this paper we do not consider such cases. 
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As follows from Eq. ([5]), the growth of the Euclidean 
norm of tangent vectors in forward-time dynamics is 
determined by the matrix ^(ii, i2)'^-^(ii, ^2)- We de- 
note its eigenvectors and eigenvalues as ff{ti,t2) and 
cri{ti,t2)'^, respectively, where cri(ii,i2) > 0'2(ii,t2) > 
• • • ^ <^m(ti,t2) > 0. The eigenvectors are termed opti- 
mal vectors because the maximal growth ratio is equal 
to (Ti(ii,t2) and is achieved if the initial vector v(ti) co- 
incides with fi{ti, 12). The same role for the backward- 
time dynamics plays the matrix t2)~'^-^(ii, ^2)"^ 
with the reciprocal eigenvalues and the eigenvectors 
f-{tiM)- 

The eigenvectors and eigenvalues can be found via sin- 
gular value decompositions (SVD) [23| of the propagator 
matrix and its inverse, thus: 

:F{hM)ft{hM) = n{tut2)<y^{tut2), (7) 
HhMr^nihM) = ^2)^1,(^1,^2)-^ (8) 

Here (Ji{ti,t2) are called singular values, and ff{ti,t2) 
and f~{ti,t2) are right and left singular vectors of 
^(ti , ^2), respectively. The singular vectors are orthonor- 
mal. They are norm-dependent, i.e., they have different 
orientations with respect to different norms 0, fill . Tak- 
ing into account Eqs. ([S]), ([7]), and ([5]), one can write the 
SVD for the adjoint propagator ^(^1,^2) and its inverse 
as 

G{tl,t2)ft{tl,t2) = f~{tl,t2)<J^{tut2)-\ (9) 

Q{tiM)~^n{tiM - ft{tiM)a^{tut2). (10) 



Comparing Eq. ^ with jS]) and Eq. ^ with ^ we see 
that the propagators JF(<i, ^2) and Q{ti, 12) have identical 
singular vectors and reciprocal singular values. 

If all ai(ti,t2) are distinct, the singular vectors are 
unique up to a simultaneous change of signs of elements 
of ff{ti,t2) and f~{ti,t2). In the presence of degener- 
acy, we still can find a set of orthonormal right singular 
vectors that are mapped according to Eq. ^ onto a set 
of orthonormal left singular vectors, but these sets are 
not unique and can be selected arbitrarily. 

Strictly speaking, propagators and singular vectors as 
well as the Lyapunov vectors considered below can de- 
pend on time both explicitly, and implicitly via state 
vectors u{t). To avoid complicated notation, we shall 
use a compact form, like .^(^1,^2)- 



B. Properties of propagators. Transformation of 
volumes built on singular vectors 

Let us discuss how ^{ti,t2) transforms volumes of 
different dimensions: segments, squares, cubes and so 
on. Being at a trajectory point at ti we construct a k- 
dimensional unit volume using the first k right singu- 
lar vectors ff{ti,t2). According to Eq. ([7]) T{ti,t2) 
transforms these vectors into the left singular vectors 
/i"(^ii^2) associated with the trajectory point at t2 



that are stretched/contracted by factors (Ti(t 1,^2), see 
Fig.[lja). The volume at t2 is equal to the product of the 
first k singular values. Alternatively, we can consider a 
fc-dimensional ball of unit radius at ti. At <2 this ball is 
transformed into an ellipsoid with axes along the vectors 
and lengths Ui. One can describe this transformation 
of volumes by 



Vk{t2) = Vk{ti)exp { {t2-ti)}2HtiM 



(11) 



where Vk{t) is the A:-dimensioiial volume, and jli{ti,t2) = 
Inai{ti,t2)/{t2 — ti) are stretch ratios that can be consid- 
ered as local Lyapunov exponents. (Note that there are 
alternative definitions of local Lyapunov exponents that 
shall be considered below.) 

The backward transformation with J^{ti, ^2)^^ is sym- 
metric. At t = t2 we construct a unit volume using 
the first k left singular vectors f^{ti,t2). According to 
Eq. ([8]) , the right singular vectors span this volume at i = 
ti, and the edges of this volume are stretched/contracted 
by factors ct~^ , see Fig.[Tfb). In a similar manner we can 
consider a unit ball at t2 that is transformed into an ellip- 
soid at ti. Therefore, the volumes are again transformed 
in accordance with Eq. PT]) . 

This discussion is also valid for the adjoint propaga- 
tor 5(^1,^2)- But because the singular values are now 
reciprocal, the volumes are transformed as 



Vk{t2)^Vkiti)cxp{ ~{t2 



tl)^fM{tl,t2) 
i=l 



(12) 



C. Far-past and far-future operators. Forward and 
backward Lyapunov vectors 

For infinitely large time intervals we can expect to ob- 
tain limits for the stretch ratios and singular vectors. 
The Oseledec multiplicative ergodic theorem [2l| and its 
corollaries state that the limit indeed exists for t2 ^ 00, 
and also a limit can be reached for ti — )- —00. When 
^2 00, the far-future operator is defined as 



W+{t)= hm [j^{t,t2fHt,h)] 



l/(2(t2-t)) 



(13) 



= lim 

io— >CXD 



F+(t,t2)S(i,t2)'/(*^-*)F+(t,t2) 



where F+(t,t2) = [/^ (t, ^2), • ■ • , /+ (t, ^2)] andS(t,i2) = 
diag[cri (t, ^2), ■ • ■ , cTrnit, ^2)] are matrices of singular vec- 
tors and values, respectively. The eigenvectors of the 
far-future operator are the limits of vectors f^{t,t2). 
We denote them as ip^{t) and refer to them as forward 
Lyapunov vectors. They are orthonormal and depend on 
t [7| ■ The convergence of the singular vectors to the Lya- 
punov vectors is considered in Ref. [l^. Logarithms of 
eigenvalues of W^(t), Ai > A2 > • • • > Xm, are called 
Lyapunov exponents. Regardless of time dependence of 
{t) , they do not depend on time. 
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Figure 1. Transformation of a volume, (a) Forward step by the propagator J^{ti,t2), (b) backward step via J-{ti,t2) 



The far-past operator is defined as 



ti—^ — oo 

= lim 

ti^ — oo 



(14) 



where F"(ii,t) = [/^(ti, i), ■ • • , /m(ii, 0]- The eigen- 
vectors of this matrix are the Hmits of the left singular 
vectors f^{ti,t) for ti —oo. They are called back- 
ward Lyapunov vectors. These vectors are also referred 
to as Gram-Schmidt vectors, because they can be com- 
puted in the course of a procedure, which includes Gram- 
Schmidt orthogonalizations; see. Sec.|lll We denote them 
by ip^{t). Similar to the forward vectors, the backward 
Lyapunov vectors are orthonormal, and depend on < 0|. 
As well as singular vectors, forward and backward Lya- 
punov vectors are norm-dependent 0, [TT| . The loga- 
rithms of the eigenvalues of W~ (t) are equal to the Lya- 
punov exponents with opposite signs. 

In analogy with the finite-time case, the /c-dimensional 
volumes can be built on the forward Lyapunov vectors 



ip^{t). Modifying Eq. (fTT|) we find that average growth 
rates of these volumes are the sums of Lyapunov expo- 
nents. 



fe 

E 



Xi = lim 

t2— >00 



1 



t2-t, 



■In 



Vk{t2) 

Vk{ti) 



(15) 



As we shall see below, this formula is valid for almost 
any fc-dimensional volume in the tangent space, not nec- 
essarily related to the forward Lyapunov vectors. 

The Lyapunov exponents may not be all distinct. 
To take possible degeneracy into account we introduce 
an additional notation. Let s be a number of dis- 
tinct Lyapunov exponents {I < s < m), and let A*-*' 
(i = 1, 2, . . . , s) denote the ith distinct Lyapunov expo- 
nent with the multiplicity So, we have A*-^^ > A*-^^ > 
• • • > A*^''-* , and 'J2i=i — what follows, to ad- 

dress the whole spectrum of Lyapunov exponents as well 
as related vectors, we shall employ lower indices while 
paying special attention to the multiplicity, we shall use 
upper indices. The notation ^^(i-, will stand for a set of 
vectors, related to the ith distinct Lyapunov exponent. 
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and where j = 1,2, ... , v^^\ will denote the jth 

vector related to A^*^. 

In presence of the degeneracy forward and backward 
Lyapunov vectors are not unique. But as we already 
mentioned for singular vectors, this is not an obstacle, 
because it is always sufficient to choose any orthonormal 
set of these vectors. 

The adjoint propagator Q can also be used to define 
far-past and far-future operators and forward and back- 
ward vectors, respectively. The Lyapunov exponents are 
the logarithms of the eigenvalues of the far-past opera- 
tor, while the far-future operator is associated with the 
Lyapunov exponents with inverted signs. 



D. Oseledec subspaces. Asymptotic behavior of 
arbitrary vectors and volumes 

Let us now discuss what happens with arbitrary vec- 
tors. The framework that helps to understand it is pro- 
vided by the following set of subspaces: 

S+{t) ^ span = J, J + 1, ■ • • , s} , 5+ ^(i) = 0, 

St (t) C SUit) c • • • C 5+ (t) = R"\ (16) 

In other words, S^it) is spanned by forward Lyapunov 
vectors {t) [i > j) related to the distinct Lyapunov 
exponents starting from the jth one. Dimensions of these 
subspaces are dim Sj'{t) = Xli=j ; where j/*^*^ is the 
multiplicity of A'-'' . Analogous subspaces spanned by the 
backward Lyapunov vectors (t) are defined by 

Sr{t) = span {cp-,., (t)|i = 1, 2, . . . ,j} , S^{t) = ill, 

S^{t)cS^{t)c---cS;{t)^W"\ (17) 

and their dimensions are dim S-{t) = T.Ui"'-'^- These 
sets of subspaces are referred to as Oseledec splitting 0, 

MM- 

Recall that the propagator lF{ti,t2) maps each right 
singular vector onto the corresponding left singular vec- 
tor and stretching rates are determined by singular val- 
ues, see Eq. ([7]). When (^2 — ^i) — >■ oo, the right and left 
singular vectors converge to forward and backward Lya- 
punov vectors, respectively, and the stretching rates con- 
verge to the Lyapunov exponents. Hence, the Oseledec 
subspace [t) consists of vectors that asymptotically 

grow or decay with rate A < A'-'^ . In turn, the vectors 
from Oseledec subspace Sj (t) grow or decay with expo- 
nential rates A > A'--'^ backward in time. 

Consider a vector t;(-')(i) G 5+(t)\S'+|.i(t). This vector 

is orthogonal to each tp^ii) (^)' where i < j, and obligatory 
has a nonzero projection onto at least one of the vectors 
ftij) (t), related to the jth distinct Lyapunov exponent. 
It means that being iterated for infinitely long time with 
the propagator the vector v''^\t) exponentially grows 



with the average rate A^^^ (Ml, [H, [H] , 

\\T{h,h+t)v^^\h)\\^c^'''K (18) 

The vectors (t) G Sj{t) \ Sj_-j^{t) behave analogously 
in backward time: 

\\J^ih-t,h)-'v^^Hh)\\^e-^'''K (19) 

Vectors v^'^'>(t) g S^it) \ S^{t) fill almost the whole tan- 
gent space, because the excluded subspace S^it) has a 
measure zero in M™ . It means that under the action of 
almost any vector, i.e., 1-dimensional volume, asymptot- 
ically grows or decays with the exponent A'^' , and its im- 
age tends to the subspace span{y>~(i, (t)} = S'f (t). Con- 
sider now a square, i.e., a 2-dimensional volume. First 
we assume that A*-^^ is not degenerate so that i/^^^ = 1. 
Almost any such square has a 1-dimensional intersection 
with the subspace S^it) \ Stit) of vectors v^^\t) that 
are dominated by the A^^) 0, [H-ill. (Here "almost" 
means that there is a measure zero set of squares fully 
belonging to subspaces with j > 1.) Thus, the area of the 
square asymptotically grows or decays with the exponent 
A'^-'^^-l-A^^^ . All segments within this square except a single 
one approach the subspace span{y~(i) (t)}, while that one 
goes into span{y~(2) (i)}. As a result, this square tends 
into the subspace 5^. When z^'^' = 2, the area of the 
square grows/decays with 2A'^'=Ai-|-A2 and the whole 
square is embedded into . But when we take a cube, 
its volume grows or decays with 2A(i)+A(2) = A1+A2 + A3 
and its image goes into S2 ■ In general this can be 
formulated as follows. Under the action of ^ almost 
any fc-dimensional volume asymptotically grows or de- 
cays with average exponential rate X^iLi ^^'^ tends to 
settle down inside the subspace S~ , where i is defined 
from the inequalities dimS'^^j^ < k < dim 5". In the 
same way considering vectors t;^-'^(<) £ S~ {t)\S~_i{t) we 
see that almost any /c-dimensional volume being iterated 
in backward time with the propagator J^{ti,t2)~^ grows 
or decays with the exponential rate X]i=i ^-'^d settles 
down in Sf{t), such that dimS'+_;^(t) < A: < dimS'+(i). 
Formally, these asymptotic embeddings can be described 
as: 

J^{h,t)Vk[h) C Sj[t), 

J ^ ^20) 
dim Sj_i{t) < k < dimS'^^(t), 

:F{t,t2)-'Vk{t2) c sUt), 

dim S'+^i(t) < fc < dimS'+(t). 

Let us now turn to the adjoint propagator G(ti, ^2)- We 
recall that its singular vectors coincide with the singular 
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vectors for J- ^ while its singular values are reciprocal. 
Hence the adjoint Oseledec subspaces can be defined as 

B^it) = span (t) |* - 1, 2, . . . , j} , (t) 0, 

Bt(t) (lH+{t)(l---Cl H+{t) = R", (22) 

Hr{t) = span {t)\i ^ j,j + 1, . . . , s} , H;^,{t) ^ (d, 

H;{t) C H;_,{t) C • • • C i/r(i) = K'"- (23) 

Note that H+_j^{t) ± 3^(1) and i?7+i(i) -L 5'-(t). Rea- 
soning in the same way as above, we find that the adjoint 
propagator Q generates the following asymptotic behav- 
ior as t oo: 

v^^\h)eH+{U)\H+_,{t,)^ 

\\gih,h+t)v^'\h)\\^e-^'''', (24) 

v<~^\h)eH-{t,)\Hr^,{t,)^ 

\\g{h-t,hr'v^^\t,)\\r^e^'''\ (25) 

and the asymptotic embeddings read: 

giti,t)Vk{h) c H^{t), 

ti^-oo (26) 
dirnHj-^^it) <k< dimiJr(t), 



git,t2)-'v„it2) c H+{t), 

dimi/+_i(t) < A; < diniiJ+(i). 

E. Finite-time evolution of forward and backward 
Lyapunov vectors 

Now we need to discuss how orthogonal Lyapunov vec- 
tors are transformed in finite time intervals. First con- 
sider the action of ^(^1,^2) on forward Lyapunov vec- 
tors. For any such vector related to the jth distinct Lya- 
punov exponent A*--*^ we can write: -(ti) & Sj'{ti) \ 

Sj'^iiti), where i = 1, 2, . . . , see Eq. (|T6|) . It means 
that this vector shows the asymptotic behavior (|18p . i.e., 
it grows or decays with the exponent A'-'^ forward in time. 
In turn, it means that T{ti,t2)<p'^(j) -(ii) = ^(^2) G 
S^{t2) \ 5++i(i2)- We see that the image of ip^^^^ .{h) 

at t2 is orthogonal to vectors ¥'^(„)(t2) with n < j. But 
this is not a forward Lyapunov vector anymore, because 
the subspaces span{y>^(^) (t2)} and Sj'{t2) \ 5'^^(i2) are 
not identical. Vectors from S^{t2) \ S^_^-^{t2) obligatory 
have a nonzero projection inside span{(p^(j) (^2)} but typ- 
ically do not belong to it and also have projections onto 
forward vectors with n > j. 



Let us first assume that there is no degeneracy i.e., 
all Lyapunov exponents are distinct. In matrix form we 
have :F(ii,t2)*+(ii) = *'*'(i2), where 

*+(t) = [v+(t),cp+(0,...,<p+(i)] (28) 

is the matrix consisting of the forward Lyapunov vectors. 
According to the above discussion of '4'~\u) ^(^2), the first 
vector-column of ^^(^2) is coUinear with (^2)- The 
second one is orthogonal to ip^{t2), but can have nonzero 
projections onto all others forward vectors. The third one 
is orthogonal both to ip^{t2) and to y>J(t2) and so on. 
Thus we can write 

*+(t2) = *+(i2)L, (29) 

where L is a lower triangular matrix. 

When the spectrum of Lyapunov exponents is degen- 
erate, the matrix ^'^{12) is not unique. There exist 
subspaces span{c^^(j) (^2)} corresponding to each unique 
Lyapunov exponent, such that any vector from these sub- 
spaces can be treated as a forward Lyapunov vector. This 
means that the decomposition (|29p is also not unique, 
because there exists a variety of non-triangular matri- 
ces L satisfying this equation. But the representation of 
'^'^(^2) as a product of an orthogonal and a lower tri- 
angular matrices exists and is unique regardless of the 
degeneracy of Lyapunov exponents. In fact, this is the 
well-known QL factorization The analysis of the 

details of the factorization procedure shows that the or- 
thogonal matrix can always be treated as a matrix of 
forward Lyapunov vectors. Hence, regardless of the de- 
generacy, Eq. ([29|) remains valid. 

Altogether, the propagator ^ maps forward Lyapunov 
vectors onto new vectors that are not Lyapunov vec- 
tors. In other words, forward Lyapunov vectors are non- 
covariant with the dynamics. To recover forward Lya- 
punov vectors, we have to perform a QL factorization. 
For the subsequent analysis it is convenient to represent 
it as a mapping backward in time: 

:F{h,t2r'^+it2) = *+(tl)L^(il,t2), (30) 

where L'^(ti,t2) G jj^x"' ig a lower triangular matrix. 
Because the propagator is non-singular and QL factoriza- 
tion is unique (if one requires for all diagonal elements 
of L"^(ii,i2) to be positive), this equation determines 
*+(ti) via *+(t2) in a unique way. By definition, the 
diagonal elements of L'^(ti,t2) do not vanish, i.e, this 
matrix is non-singular. 

Repeating the above discussion for the backward Lya- 
punov vectors, we see that regardless of the degeneracy, 
the following relation is always valid: 

J^{h,t2)^-ih) = *-(t2)R^(ti,i2), (31) 

where 

^-{t)^[v^{t),^^{t),...,^-M, (32) 
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Figure 2. The idea of orthogonalization. The vectors Vi are the result of mapping ([5]) and vectors q- are their orthogonalization: 
qj is collinear to vi, qj belongs to the plane spanned by vi and V2, and q,^ belongs to the space spanned by vectors vi, V2, 
and V3. 



and R"^(ti,i2) is an upper triangular matrix with a 
nonzero diagonal. Like the forward vectors, the back- 
ward vectors are non-covariant with the dynamics. 
For the adjoint propagator ^(^1,^2) we obtain: 

G{h,t2)-'^+{t2) = *+(ti)R^(ii,i2), (33) 
g{h,t2)^-{h) = *-(t2)L^(ti,t2), (34) 

where R^(ti,t2) and Li^{ti,t2) are upper and lower non- 
singular triangular matrices, respectively. 

II. NUMERICAL COMPUTATION OF 
LYAPUNOV EXPONENTS AND FORWARD AND 
BACKWARD VECTORS 

The definition of the Lyapunov exponents and vectors 
cannot be implemented directly as a numerical algorithm. 
It is impossible to solve Eq. ([3]) for a sufficiently long time 
interval t2~ti, to calculate the propagator J^{ti, 12), and 
then to find a good approximation for the limit matrix 
W^. As we already discussed above, when we move away 
from the starting point ti almost any vector approaches 
the first backward Lyapunov vector (pj~(t), i.e., falls into 
subspace S^{t). Hence, in this way we can compute only 
the largest Lyapunov exponent and the corresponding 
vector. 

Equation pip determines a mapping of backward Lya- 
punov vectors at ii onto backward Lyapunov vectors at 
t2. A set of all backward vectors at different times can 
be considered as a kind of limit set, attracting or re- 
pelling, and the mapping (|3T|) can be treated as sta- 
tionary dynamics on this set. This gives an idea for 
an iterative computation of the backward Lyapunov vec- 
tors. One can initialize an arbitrary orthogonal ma- 
trix and start iterations including mapping by ^ and 
QR factorization as described by Eq. ([5T|) . These itera- 
tions converge to the backward Lyapunov vectors where 
convergence is guaranteed by Eq. (PO)) . One sees that 



the forward in time mapping embeds an arbitrary vol- 
ume into the subspace spanned by backward Lyapunov 
vectors. It means that in the course of forward itera- 
tions T{tn,tn+i)Cl{tn) = Q(t„+i)R-^(t„+i , t„) columns 
of Q{tn) € R™'^'" converge to backward Lyapunov vec- 
tors. In fact this idea was suggested almost simultane- 
ously by Benettin et al. (23l . l2a | and by Shimada and 
Nagashima [23| to compute the Lyapunov exponents. 
The convergence of these iterations towards the back- 
ward Lyapunov vectors is discussed in Refs. 0, Q- 

Consider the iterations in more detail, see Fig. [2l Sup- 
pose we have an orthogonal matrix Q(tn). First we de- 
termine !F{tmtn+i) for some interval tn+i — tn, which 
typically is not very large, and perform the mapping 
V(i„+i) — ^{tn,tn+i)Q{tn)- The first vector-column 
Vl of V(t„_(_i) behaves as we need, namely it approaches 
the subspace . So, we only normalize it to prevent 
overflow or underflow: vi qi, \\qi\\ = 1. The plane 
spanned by vectors vi and V2 approaches the subspace 
S2 if Al 7^ A2, or it goes into otherwise. In the first 
case we need to prevent the collapse of the plane due to 
the alignment of V2 along tp^ , and also the orientation 
of the plane has to be preserved to support the conver- 
gence. These two goals can be achieved by finding a new 
vector ^2 which is orthogonal to and belongs to the 
plane originally spanned by Di and V2. This vector is 
also normalized. In the second case, when Ai = A2, there 
is no alignment and, in principle, there are more options 
how to define But it is allowed anyway to compute 
^2 as if the degeneracy was absent, and this is the most 
reasonable choice making the procedure most transpar- 
ent. In a similar manner we find the third normalized 
vector that is orthogonal to qi and ^2 and belongs to 
the space spanned by Vi, V2 and D3. Doing so for all the 
remaining columns of V(f„-|_i) we compose the matrix 
Q{tn+i) whose columns are vectors q^. Then we use this 
Q(f„+i) as an initial value for the next mapping with 
^(tn+i,tn+2) and repeat the procedure. After many re- 
cursions the columns of Q{tn) converge to the backward 
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^{tn, tn + l) 




ri2 Qi vi = 



Figure 3. Computation of a volume after the mapping by J-{tn, tn+i)- 



Lyapunov vectors. This procedure works not only for the 
whole set of vectors, but allows to compute any number 
of the first backward Lyapunov vectors. 

The described procedure eliminates the ambiguity of 
backward Lyapunov vectors that emerge when not all 
Lyapunov exponents are distinct. Particular directions 
of backward Lyapunov vectors corresponding to each de- 
generate Lyapunov exponent A'-'' depend on the choice 
of the initial matrix Q(to)- But these variations remain 
within subspace span{y } so that any choice is appro- 
priate. Moreover, in practical computations the degen- 
eracy manifests itself very weakly, because typically the 
degenerate Lyapunov exponents converge to identical val- 
ues very slowly. In fact, dealing with a high-dimensional 
system one needs to know in advance which of the expo- 
nents are expected to be identical to identify them in the 
computed spectrum. 

The computation of the Lyapunov exponents is illus- 
trated in Fig. [21 An initial unit square composed of vec- 
tors Qiitn) is transformed into the parallelogram spanned 
by the vectors Vi{tn+i). After the orthogonalization we 
obtain q^{tn+l)■ To compute the area of the parallel- 
ogram we can construct a rectangle with identical area 
by projecting Vj{tn+i) onto qj(<„+i): = q^v^. As 
we see from the figure, the area is riir22- Similarly, a 
A:-dimensional unit volume after the mapping is equal to 
■ ■ - fkh- Thus, we can define the local Lyapunov 
exponents as 

Xi = In(rij)/ {tn+i - t„). (35) 

In the course of the mapping/orthogonalization iterations 
we need to accumulate and average Ai to obtain the Lya- 
punov exponents. 

By construction, the first vector Vi has only one 
nonzero projection onto Qi, the second vector V2 has 
two nonzero projections, onto Qi and the third vec- 
tor i>3 has three nonzero projections onto first three 
vectors and so on. It means that Vij are elements 
of an upper triangular matrix. So, the procedure de- 
scribed above represents the matrix V as the product 



V = QR. Here Q is an orthogonal matrix such that 
span{qi,q2, ■■■Qk} ^ span{t>i, 1)2, . . .ffc} for any k < m, 
and R is an upper triangular matrix consisting of the pro- 
jections of columns of V onto columns of Q. This proce- 
dure is called QR factorization There are different 
numerical algorithms of the QR factorization. Note that 
the often used Gram-Schmidt algorithm as well as its 
modified version are not very accurate when the dimen- 
sion of the tangent space is large Most high preci- 
sion QR algorithms are based on so called Householder 
transformations [20I [28j. 

Another way to compute backward Lyapunov vectors 
is based on the adjoint propagator Q. Equation ([34]) de- 
termines the stationary dynamics, and Eq. ([^H) indicates 
that the forward iterations converge to this dynamics. 
Because Q{tn,tn+i) has reciprocal singular values, the 
value cr„j(i„, t„_|-i)~^ dominates in the course of forward 
iterations with the adjoint propagator. It means that 
columns of Q converge to the backward Lyapunov vec- 
tors in the reverse order. If we rearrange columns of $~ 
in Eq. ([34| in the reverse order, we also have to trans- 
pose with respect to its diagonal and with respect 
to the antidiagonal. As a result we obtain an upper tri- 
angular matrix. Thus, the algorithm is identical to the 
one previously described. We perform the mapping by 
Q{tn,tn+i), find a QR factorization of the resulting ma- 
trix, take Q{tn+i), and do the next recursion. 

Consider now the computation of the forward Lya- 
punov vectors. The first algorithm is based on Eqs. ((?7|) 
and ([33)) . We need to move backward in time alternat- 
ing mappings with Q{tmtn+i)~^ and QR factorizations. 
The matrices Q converge to and the forward Lya- 
punov vectors come up in the correct order. Note that 
is merely the transposition of JF, see Eq. ([6|). In 
the course of this procedure we can compute local Lya- 
punov exponents as logarithms of diagonal elements of 
triangular matrices per unit time. For short time inter- 
vals these local exponents will differ from those given by 
Eq. (|35p , but being averaged over many times steps they 
also converge to the Lyapunov exponents. 
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Another algorithm for the forward Lyapunov vectors is 
based on Eqs. (|30p and ((2T|) . The procedure is the same 
as above except using the inverted propagator J-~ ^ . This 
method computes the vectors in the reversed order, and, 
hence, the previous one is usually more applicable. The 
idea to apply the transposed propagator instead of the 
inverted one was suggested in Ref. 

The implementation of the algorithm with the trans- 
posed propagator !F"^ is straightforward for discrete time 
systems (e.g. coupled map lattices), where the action 
of on a set of (Lyapunov) vectors can be computed 
using the transposed Jacobian matrix of the system. In 
principle, one can do the same with continuous systems, 
but in that case one would have to compute the full 
propagator first by solving m copies of the linearized 
ODE ^ and then use its transpose J^"^ to evolve the 
desired number of tangent vectors. This implementa- 
tion is inefficient if the system is very high dimensional 
(to 3> 1) and if only a few Lyapunov vectors are to be 
computed. As an alternative, the action of can re- 
formulated as follows. Using the Magnus expansion (29| . 
we can represent the propagator of Eq. ([2]) via matrix ex- 
ponential functions as jF(<i,t2) = exp[r2"^(ti, ^2)]- Here 
is a matrix that is given as a series expan- 
sion n^{h,t2) = E^i^^f (^i'^2), with nf{ti,t2) = 

J^^3in)dn, n^{h,t2) = i/t?rfn/,7dT2[J(ri),J(r2)], 
and so on, see [l^, J(t) = 3{u,t) is the Jacobian 
matrix, and [•,•] denotes the matrix commutator. The 
adjoint propagator reads: Q{ti,t2) = J^{ti,t2)~'^ ~ 
exp{-[r2-^(ii,i2)]'^} = exp[rj^(ii,i2)]- The matrix 
^^{ti,t2) = ~[n^{ti,t2)r generating ^(^1,^2) is ob- 
tained with a Magnus expansion where the Jacobian ma- 
trix 3{u,t) is replaced by —J{u,t)'^. So, to compute the 
action of ^(^1,^2) on a tangent vector we have to solve 
the following linear ODE 



V = — J(m, t)'^v 



(36) 



forward in time (from ti to t2 > ti , because the action of 
the adjoint propagator Q(ti,t2) corresponds to moving 
forward in time). To compute forward Lypaunov vectors 
using J^{ti,t2)'^ — Q{ti,t2)~^ we have to invert Q{ti,t2)- 
This can be done by integrating the required number of 
copies of Eq. ([36|) and the basic system ([ij backward in 
time (from t2 to ti). 

All four algorithms compute the dominating Lyapunov 
exponents and corresponding vectors with the highest 
precision, while the remaining part of the spectrum is 
not very accurate. Namely, IF- and ^^^-algorithms do 
the best for the first Lyapunov exponent and vectors, 
while and ^-algorithms achieve the highest accu- 

racy for the TOth exponent and vectors. One can per- 
form and ^-algorithms in parallel, and then construct 
weighted sums of computed exponents and backward vec- 
tors to obtain the whole spectrum with very high preci- 
sion. Similarly, performing backward iterations simulta- 
neously with and one can compute the forward 
Lyapunov vectors with improved accuracy. 



III. COVARIANT LYAPUNOV VECTORS 

Orthogonal matrices computed according to QR de- 
composition preserve subspaces spanned by each first k 
columns of a factorized matrix. The QL decomposition 
preserves subspaces spanned by each last k columns of 
a factorized matrix. It means that considering Eqs. ()30p 
and (|3T|) . we can conclude, that Oseledec subspaces (fT6|) 
and pT|) are preserved under the tangent flow P, 0| . The 
same conclusion follows from Eqs. ([55]) and for the 
subspaces (|22p and 



(37) 



T{tl,t2)S+{ti)=S+{t2), 

T{ti,t2)S-{t,)^Sr{t2), 

Gitl,t2)H+{ti)=H+{t2), 

g{tiM)H;{ti)=HT{t2). 



(38) 



So, the Oseledec subspaces are invariant under time re- 
versal and covariant with the dynamics. But this is not 
the case for the forward and backward Lyapunov vectors 
themselves. Being multiplied by ^ and Q they also have 
to be multiplied by lower or upper triangular matrices 
to be mapped to new forward and backward Lyapunov 
vectors, see Eqs. dSO]), jSIl), ([33]), and ((34)) . 

Given the covariant subspaces, it is natural to search 
for some vectors inside these subspaces that are also co- 
variant with the dynamics and are invariant with respect 
to time reversal. These vectors are referred to as covari- 
ant Lyapunov vectors (T^ . We denote them by 'Jj{t). 
The basic property of these vectors (which are covariant 
with respect to the propagator ^) can be written as 



\T{ti,t,±t)jJt,)\\ ^cxp(±A,t) 



(39) 



for any ti and t — >■ 00. The covariant Lyapunov vec- 
tors are norm-independent 0, [U . Also we can introduce 
norm-independent adjoint vectors 6j(t) that are covari- 
ant with respect to the adjoint dynamics: 



\G{ti,ti±t)e,iti)\\ -exp(TAji). 



(40) 



Equation ([391) means that Eqs. (HH) and ([19]) are fulfilled 
simultaneously, and Eq. (PO)) implies the simultaneous 
validity of Eqs. (|24p and ([25]) . It means that the covari- 
ant Lyapunov vectors belong to the intersection of the 
Oseledec subspaces P, 0, f30l |. and the adjoint covariant 
vectors can be found within the intersections of the ad- 
joint subspaces: 



f^{t)es+{t)ns-{t), 
e,{t)eH+{t)nH-{t). 



(41) 
(42) 



These intersections are always nonempty because the 
sum of dimensions of Oseledec subspaces is always higher 
than the dimension of the whole tangent space. 

Consider arbitrary vectors w'^-'^(ti) € Sj'{ti)\Sj'^-^{ti), 
where j ~ 1, 2, . . . , s, and s is the number of distinct Lya- 
punov exponents. There are i/^^' linearly independent 
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vectors corresponding to the jth Lyapunov exponent A^-'^ , 
and the total number of such vectors is X]j=i ~ 
Representing the whole set of these vectors as a matrix 
V, we obtain V = $~''A~'', where A"*" is a lower trian- 
gular matrix, and is a matrix of forward Lyapunov 
vectors (P5|) . As follows from Eq. (fT5|) . when the for- 
ward propagator ^ is applied to these vectors, the first 
z/^^) of them grow or decay asymptotically with the ex- 
ponent A^^-', the next v'^'^^ vectors grow / decay with the 
exponent A^^^ and so on. In a similar manner we can 
consider arbitrary vectors v^3)[t^ g '^^ \ 
The matrix of these vectors V can be found as V = 
#^A~, where A~ is an upper triangular matrix, and 
^ is defined by Eq. According to Eq. HH), act- 

ing upon these vectors by the inverted propagator ^ 
we can observe that the first j/'^^ of them grow or de- 
cay asymptotically with the exponent — A'^\ the next 
i^'^^ vectors grow / decay with the exponent — A^'^-' and 
so on. Let T(t) = [71(^)172(^)7 ••• 77m] be a matrix 
consisting of the covariant Lyapunov vectors, and let 
0(t) — \Q\(t\ Q-iifyi ■ ■ ■ 1 ^m] be a matrix of adjoint co- 
variant vectors. As follows from Eq. ([M)) . the covariant 
vectors have to demonstrate both forward and back- 
ward p9p asymptotic behavior. It means that there exist 
an upper triangular matrix A^ and a lower triangular 
matrix A^, such that 

r(t) =*-(t)A-(t) = *+(t)A+(i). (43) 

Reasoning in a similar manner one obtains for the adjoint 
vectors: 

0(t) = $+(f)B+(t) = *~(i)B~(i), (44) 

where B^(t) and B^(t) are upper and lower triangular 
matrices, respectively. Note that Eqs. (HS)) and (^1)1 con- 
vey, in fact, the same property of covariant vectors as 
Eq. (jH]) and ([42]), respectively. Muhiplying Eq. (gSl) by 
[*+(t)]'^ and Eq. dH]) by [*"(i)]'^ we obtain the rela- 
tions between triangular matrices that will be required 
later: 

P(t)A-(t)=A+(0, (45) 
P(t)TB+(t)=B-(i), (46) 

where 

P(t) = [*+(t)]'^*-(t) (47) 

is a m X m orthogonal matrix. 

If the Lyapunov exponents are degenerate, the covari- 
ant vectors are not unique. Let us discuss what Eq. (j43[) 
implies in this case (Eq. (|44p can be considered in the 
same way). If r(t) is known, then we can compute ^~(i) 
and A~(<), and *^(i) and P^{t) via QR and QL de- 
compositions, respectively, in a unique way. However, 
Eq. ([33]) does not determine V(t) via *"'"(<) and *~(i) 
in a unique way. In principle, there exist orthogonal ma- 
trices #~ and that allow one to fulfill Eq. (^5)) with 



several couples A~(i) and A~''(t), resulting in different 
matrices F, and, hence, in different covariant Lyapunov 
vectors. As an example one can consider a matrix 
that consists of columns of arranged in the reverse 
order. Ambiguity of A^(t) and A^(i) means that there 
are Lyapunov exponents associated with several covari- 
ant vectors. But on the other hand, the total number 
of covariant vectors is equal to the total number of Lya- 
punov exponents to, and there are no exponents without 
vectors. It means that the ambiguity can occur if and 
only if the Lyapunov exponents are degenerate. The co- 
variant vectors associated with a fc times degenerate Lya- 
punov exponent can have arbitrary orientation within a 
fc-dimensional subspace corresponding to this exponent. 
But because any set of linearly independent covariant 
vectors from the subspace corresponding to the degen- 
erate exponent is as good as any other, this ambiguity 
can be ignored: we just need to have any linear inde- 
pendent set of vectors. (We recall that though forward 
and backward vectors are also subject to the degeneracy, 
their ambiguity is eliminated in the course of the compu- 
tations, see Sec. HIl) 

Let us find how r(ii) is transformed by .^(^1,^2)- In 
general we can write 

:p(ti,t2)r(ti) = r(t2)c^(ii,i2)7 (48) 

where C'^(ti,t2) is a matrix whose structure should be 
determined. When the Lyapunov spectrum is not de- 
generate, Eq. (|1T|) immediately implies that C"^(ti,t2) 
is diagonal. To show that this is the case regardless of 
the degeneracy, we substitute V{\) ~ $^(t)A^(t), see 
Eq. (gSl), in Eq. (gS]) and, taking into account Eq. ([gP)) . 
we obtain 

L-^(ii,t2)A+(i2)C^(ti,t2) = A+(ti). (49) 

Since all known matrices here are lower triangular, 
C"^(ti,t2) is also lower triangular. Analogously substi- 
tuting r(t) = *"(t)A~(t) in Eq. ^ and using Eq. ^ 
we obtain: 

R-^(ti,t2)A-(ti) = A-(t2)C^(ii,i2), (50) 

i.e., C'^(ii,i2) is an upper triangular matrix. Si- 
multaneous upper and lower triangular struc- 
ture has only a diagonal matrix: C"^(ti,t2) = 
diag[ci(ii,t2),C2(<i,i2), ■ ■ • ,Cm(<i,<2)]. Hence, the 
vectors 7^ can freely evolve under the tangent flow 
so that the tangent flow preserves their directions. The 
direction, represented by 7j(ii) at t\ is mapped onto 
the direction pointed by at t2, and the backward 

step maps Ijit-i) onto the direction of ')^{t\). The 
vectors themselves are stretched or contracted by factors 
cf{ti,t2)- (Recall, that the directions of the forward 
and the backward Lyapunov vectors are not preserved). 
The adjoint vectors freely evolve under the tangent flow 
generated by the adjoint propagator: 

Q{tl,t2)&(tl) - e(t2)[C^(tl,t2)]-l. (51) 
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One can say that the vectors 7^ are covariant with the 
tangent dynamics generated by ^ and the adjoint vectors 
are covariant with the tangent dynamics of Q. This is 
the reason why these vectors are referred to as covariant 
vectors. 

Since the covariant vectors are defined up to an ar- 
bitrary length, the diagonal elements of C"^ can be de- 
fined in various ways. In particular, to fulfill Eq. ((5^ 
we should not normalize the vectors, and C"^ = I in 
this case. However, in the course of numerical com- 
putations we need to avoid overflows and underflows. 
Hence, constant lengths of 7j(t) have to preserved with 
respect to the chosen norm. In this case c^{t\,t2) = 
Il^(ii,i2)7,(ii)ll/ll7,(ii)ll, and 

In[cf(ti,t2)]/(t2-<i) (52) 

can be treated as local Lyapunov exponent. The values of 
these local Lyapunov exponents depend on the norm, but 
being averaged over many / long time intervals (ti,t2)i 
regardless of the norm they converge to the Lyapunov ex- 
ponents Aj. Consider an important particular case. As 
follows from the discussions in Sec.[lll one can build unit 
volumes using the covariant Lyapunov vectors when the 
diagonal elements of the upper triangular matrix A~ are 
equal to 1 , see Eq. (j43| . Equation ([50|) describes the dy- 
namics of A~ corresponding to the tangent dynamics of 
the covariant Lyapunov vectors. When two upper trian- 
gular matrices are multiplied, the resulting matrix is also 
upper triangular and its diagonal elements are the prod- 
ucts of the diagonal elements of the multipliers. Thus, if 
the covariant Lyapunov vectors are rescaled to preserve 
ones on the diagonal of A~. then the are equal to 
the diagonal elements of R'^, and the local Lyapunov 
exponents ((52|) coincide with those defined by Eq. ([35| : 

ln[cf(ii,<2)]/(t2-ii) =A,(ti,t2)- 

Let us now discuss what it means if covariant vectors 
merge. The phase space of dynamical systems can con- 
tain structures called "wild hyperbolic sets" that are re- 
sponsible for the existence of structurally stable and un- 
avoidable homoclinic tangencies between stable and un- 
stable manifolds. In turn, the presence of these tangen- 
cies results in formation of non-hyperbolic chaotic attrac- 
tors (sH . Since covariant vectors are associated with in- 
variant manifolds of trajectories, in points of tangencies 
the corresponding vectors become collinear [l], [ll|, [T^) [3 ■ 
The same happens with the corresponding adjoint co- 
variant vectors. Collinear vectors result in a singularity 
of the matrices r(i), and 0(t). The triangular matrices 
A^(t) and B^(t) also become singular. Note, that this 
property is time-invariant: as follows from Eq. (|48p and 
(|5ip if some of covariant vectors are identical at t = ti, 
they remain identical for all time. In practice, select- 
ing an arbitrary trajectory we almost never hit exactly 
the trajectory with the tangencies. But if a trajectory 
with tangencies exists, the arbitrarily selected orbit will 
pass infinitely close to it and we will encounter with a 
nonzero frequency ill-conditioned matrices of covariant 



vectors. Note that this is not the case for orthogonal 
forward and backward vectors, which are not affected by 
tangencies. 

Now we consider how covariant and adjoint covari- 
ant vectors are related to each other. First of all no- 
tice that given r(t), one can always compute and 
A^(t) as its QR decomposition and ^^(t) and A~''(t) 
as a QL decomposition. Then one can construct the ma- 
trix P(i) = and compute 0(t) via the LU 
method as described below in Sec. IIVBI It means that 
these two sets of vectors are not independent from each 
other. However, the mutual orientation of these vectors 
can help to recover some new data. 

Transposing Eq. ([44]) and multiplying it with Eq. (j43| , 
we obtain: B+(t)TA+(i) = B"(t)TA"(t). The left hand 
side of this equation is a lower triangular matrix, while 
the matrix on the right hand side is upper triangular. 
Hence, 

B±(t)TA±(t) = A±(i)TB±(i) = D(i), (53) 

where D(t) is a diagonal matrix. Again take into account 
Eqs. dH]) and (gS]) to write: 

0(i)^r(t) = r(t)T0(i) = D(t). (54) 

The diagonal structure of D indicates that each adjoint 
covariant vector 0^ (f), j = 1, 2, . . . , m is always orthogo- 
nal to the covariant vectors 7j(t), where i ^ j- In pres- 
ence of the tangency 7j(t) = 7j_|_i(t) the jth and the 
(j + l)th diagonal elements of D vanish, i.e., correspond- 
ing adjoint and original vectors also become orthogonal: 
"ijj_i{t) _L 9j+i{t), where z = 0, 1. It means that given 
the vectors "fiit), one can find the adjoint vectors 9j{t) 
as null vectors of the matrix consisting of all 7^ (t) except 
the jth one. Notice that even if a tangency occurs, one 
still can compute 9j{t) in this way. To find how D(t) is 
varying in time, we transpose Eq. (|48p . multiply it with 
Eq. ((51|), and take into account Eq. ([6]): T{ti)^&{ti) 
C^{h,t2)T{t2)^&{t2)[C^{h,h)]-\ Hence, = 
C-^(ti,*2)D(t2)[C^(ii,i2)]"^ (recaU that afl matrices 
here are diagonal). Altogether, the elements of the diago- 
nal matrix D are cosines of angles between corresponding 
covariant and adjoint covariant vectors. Since these an- 
gles are affected by tangencies, their time averages as well 
as their temporal fiuctuations; i.e., the first and other 
moments, can be considered as characteristic numbers 
describing the structure of an attractor. The angles are 
norm-independent, because they are defined in terms of 
covariant and adjoint covariant Lyapunov vectors which 
share this property. 

If the covariant vectors are computed with a non-ideal 
accuracy, the errors will grows in course of the tangent 
dynamics. The same is the case for the adjoint covariant 
vectors. In particular, it means that if we have found 
numerically covariant vectors at ti, we cannot compute 
them at i > via Eq. ([15)) because numerical errors 
results in the divergence from the true directions. But 
nevertheless, Pazo in Ref. [32| shows that this divergence 
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is actually sufficiently slow. Hence, Eq. (|^5)l can be used 
to find an estimate for the covariant vectors at t > ti 
when t — ti is not very large. 

The covariant Lyapunov vectors are defined locally, ac- 
cording to Eqs. and and asymptotically, as 
follows from Eqs. and (1^0)) . These equations pro- 
vide two basic ideas for computing these vectors. The 
first one is to find backward and forward Lyapunov vec- 
tors for some point of the trajectory and compute an 
intersection of corresponding Oseledec subspaces. The 
straightforward implementation of this approach though 
possible, takes a lot of computational resources. We dis- 
cuss it in Sec. lTVAl In Sees. |IVB] and [IVC] more "clever" 
implementations are considered. 

The second approach is to try to arrive at asymptotic 
behavior described by Eq. ([55]) or If we initialize 

a vector, satisfying Eq. (fT9| and start iterations back- 
ward in time, after a long time we closely approach the 
limiting vectors that evolve as J^{tn,tn+i)~^Vj[tn+i) = 
Vj[tn)Cj{tn, tn+i)^^ , whcrc Cj{tn, ^n+i) are related to the 
local Lyapunov exponents (|52l) . This equation is re- 
versible, so that when the limit is reached, we can turn 
forward and arrive the opposite limit too. It means that 
the limiting vectors Vj found in this way satisfy Eq. (j39p 
and coincide with 7^. The forward iterations defined 
by Eqs. ([T8|) also converge to the covariant Lyapunov 
vectors. Similarly, the iterations initialized according to 
Eqs. ((24|) and ((25| converge to the adjoint covariant vec- 
tors. The straightforward numerical implementation of 
this approach is impossible. Due to numerical noise, vec- 
tors cannot be initialized exactly as required, and the 
numerical routines always converge to the single domi- 
nating vector. But a way to avoid this obstacle is known, 
and we consider it in Sec. IIVDI 



IV. NUMERICAL METHODS FOR 
COMPUTING COVARIANT LYAPUNOV 
VECTORS 

A. Intersection of Oseledec subspaces 

A straightforward way to find covariant Lyapunov vec- 
tors is based on Eq. (|4T|) . Given forward and backward 
Lyapunov vectors, one can construct intersections of the 
Oseledec subspaces and find the covariant vectors. To 
compute the intersection of two subspaces one can com- 
pute so called principle angles between subspaces (20l[33t . 
In brief, this method is associated with computation of 
the singular values and vectors of submatrices of the ma- 
trix (|i7)) . 

To compute the jth covariant vector one needs the first 
j backward vectors and m — j + 1 last forward vectors. 
The first backward vectors can be computed in the course 
of the iterations with the propagator and the last 
forward vectors are the result of the iterations with the 
inverted propagator !F~^ , see Sec. |lll 

Regardless of j, m+l forward and backward Lyapunov 



vectors are always required. So, this method is applicable 
for computation of the whole spectrum, but this is not an 
efficient approach if one needs only a few first covariant 
vectors. Because the forward Lyapunov vectors are com- 
puted in the reverse order, this method has a "fiattened" 
accuracy along the spectrum: the backward vectors have 
higher accuracy in first part of the spectrum, and the 
forward one are more accurate in its last part. So, the 
resulting covariant vectors have approximately the same 
accuracy for the whole spectrum. 



B. Method of LU factorization 

It is possible to avoid computation of the whole spec- 
trum of the forward or backward Lyapunov vectors to 
get only a few first covariant vectors. Two original ideas, 
which were reported in Refs. (TTI . [T2| . are discussed in 
Sees. live) and ITVDI In the current section we present a 
new approach to this problem. 

Consider Eq. ([45|) . Matrices A"*" and A~ are lower and 
upper triangular, respectively. If A~ is non-singular, we 
can rewrite Eq. (|15|) as P = A^(A~)~^. This equation 
can be considered as an LU factorization of P, i.e., rep- 
resentation of a matrix as a product of a lower and an 
upper triangular matrix [23|. If the factorization exists, 
it is unique up to the diagonal elements of one of the ma- 
trices (factors). For us it means that if we find the LU 
decomposition of P, we find the covariant vectors up to 
arbitrary lengths. 

There are many well developed standard routines com- 
puting the LU factorization. But for us the serious dis- 
advantage is that they work well only as long as the as- 
sumption of non-singularity of A~ remains valid. If ma- 
trices A^ are singular, the straightforward factorization 
of P does not exist. The standard routines for LU de- 
composition avoid this obstacle performing preliminary 
permutations of rows and columns of P. This is not suit- 
able for us, because the order of rows and columns in P 
is essential. Moreover, the standard routines find both 
A~, and A^, while it is enough for us to have only A^. 

Let us return to Eq. (j^S]). We shall demonstrate now 
that the required elements of A~ can be found from this 
equation regardless of a possible singularity of A^. To 
compute the jth covariant vector we need to find the top 
j elements of the jth column of A~. This fragment of 
the column can be denoted as A^(l : The remain- 
ing fragment A~(j + 1 : m,j) contains zeros. Note that 
here we omit the time dependence and use parentheses to 
indicate submatrices. The matrix equation for nonzero 
elements reads: P{1 : j,l : j)A~ {1 : j, j) = A+(l:j,j), 
where P(l : j, 1 : j) is the top left square submatrix of P. 
Because A^ is lower triangular, the fragment A^(l 
of its jth column contains zeros except for the diagonal 
element A'^{j,j). As already mentioned above, the LU 
decomposition is unique up to diagonal elements of one of 
the matrices. It means that we can eliminate the equa- 
tion, corresponding to the jth row of P(l:j, l:j) and 
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write the following homogeneous matrix equation 

P(l:i-l,l:j)A-(l:j, j) = 0. (55) 

This equation allows to compute nonzero elements of the 
jth column of as the null space of the rectangular 
submatrix P(l : j — 1, 1 : 7). To obtain covariant unit vec- 
tors the solutions have to be normalized. 

Equation (|55p can, in principle, have multiple solutions 
for (1 : j, j). (In this case the rank of P(l : j — 1, 1 : 
is less than {j — 1).) As we discussed above, this ambigu- 
ity can occur only due to the degeneracy of the Lyapunov 
exponents, and we can arbitrarily choose one of the mul- 
tiple solutions. 

As follows from Eq. (|46|) , the adjoint covariant vectors 
can be computed analogously, using the equation 

(pT)(l:j-l,l:j)B+(l:j,.7) = 0. (56) 

Let us now consider the submatrix P(l : j, 1 : j). If this 
is singular, then Eq. ([55]) provides for the {j+l)th column 
the solution A~(l :m,j + 1) = A~(l ■.■m,j), i.e., the jth 
and (j -|- l)th covariant vectors coincide. The inverse is 
also true, and, hence, the singularity of the submatrix 
P(l:j, l:j) is a sufficient and necessary condition for 
merging of the jth and (j -I- l)th covariant vectors. 

As discussed above, the merging of covariant Lya- 
punov vectors indicates tangencies of invariant mani- 
folds of an attractor that, in particular, occur when 
the attractor is chaotic and non- hyperbolic [3l|. To de- 
tect the violation of hyperbolicity, one usually studies a 
distribution of angles between expanding and contract- 
ing subs pac es spanned by corresponding covariant vec- 
tors [13. Il4l. [3^ . [35! . (Another method for a numerical 
test of hyperbolicity, which does not employ covariant 
vectors, is based on the so called cone criterion (36j.) 
Analyzing properties of submatrices of P one can test 
for hyperbolicity without explicit computation of covari- 
ant vectors. Let the number of positive Lyapunov ex- 
ponents be k. Moving along a trajectory, we need to 
compute some characteristic number whose small value 
indicates the nearness of P(l : fc, 1 : A:) to singularity. It 
can be, for instance, the determinant or the smallest sin- 
gular value. A small characteristic number means that 
the trajectory passes close to the tangency. So, if the dis- 
tribution of characteristic numbers computed for many 
trajectory points is well separated from the origin, then 
the chaos is hyperbolic, and if it approaches the origin 
violations of hyperbolicity occur. 

One can also study the statistics of nearness to sin- 
gularity of all submatrices P(l:j, l:j), where j = 
l,2,...,m — 1. This can provide detailed information 
concerning properties of various limit sets embedded in 
an attractor. 

Another way to characterize an attractor is to com- 
pute the matrix D containing cosines of angles between 
covariant and adjoint covariant vectors. As discussed 
above, each merged couple of vectors, i.e., each tangency, 
is represented as a couple of zeros of the corresponding 



matrix elements. To compute D, first we find the ma- 
trix A^, then using Eq. (|45l) compute only the diago- 
nal elements of A"*", and after that compute B"*" using 
Eq. (Though only its diagonal elements are re- 

quired, we cannot get them without computing the rests 
of the columns.) Finally, we obtain the elements of D 
as products of diagonal elements of A^ and B"*"; see 
Eq. Note, that it is not required to compute the 

whole matrix D. The method allows one to find only a 
few first elements. 

Normally, one has to compute the covariant Lyapunov 
vectors for a series of subsequent points of a trajectory. A 
practical implementation of the algorithm in this case can 
be the following. We start the procedure for Lyapunov 
exponents forward in time including the iterations with 
^{ti, and QR factorizations, and perform it as long as 
required for the orthogonal matrices Q(t) to converge to 
the matrices of the backward Lyapunov vectors $~(t). 
Denote the end of the preliminary stage as ^a- After 
this point the iterations are continued, but now we store 
trajectory points of the basic system and the backward 
vectors #^(t„), see the diagram in Fig. |4l[a). The dura- 
tion of this stage depends on the number of points where 
we need to know the covariant vectors. At we stop the 
storing of $^(t„) and, moreover, stop the procedure for 
Lyapunov exponents and continue to solve only the basic 
system saving the trajectory points. This stage lasts from 
te to tc- Its duration must be long enough for the sub- 
sequent backward procedure to converge. At tc we start 
moving back along the saved trajectory performing the 
backward procedure for Lyapunov exponents including 
iterations with the adjoint propagator and QR fac- 
torizations. Upon the arrival at we have the forward 
Lyapunov vectors 4>^(t). Now we pass the interval from 
ts to given both the backward vectors $"(<„), that 
were saved in the course of the forward pass, and the 
forward vectors ^~^{tn). These vectors can be used to 
compute the matrices A~(t„) by means of P (Eq. ((47l) ). 
as explained before. In turn, these matrices can be used 
to find the covariant vectors r(t„), according Eq. (|43| . 
Note that it is not necessary to perform this procedure 
with the whole set of vectors. To compute j first covari- 
ant vectors we need j first backward vectors and j — 1 
first forward vectors. In the appendix we provide a pseu- 
docode implementation of the presented algorithm. 

Columns of A~(t„) can also be considered as covari- 
ant Lyapunov vectors written with respect to the basis 
^~{tn). The covariant vectors in the form of A~(i„) 
have a mutual orientation that is identical to T(tn)- 
Therefore, if for example the angles between covariant 
Lyapunov vectors are required, they can be computed 
with respect to columns of A~(t„). This allows us to 
save some machine time. 

The numerical implementation of the described pro- 
cedure includes well established numerical routines. To 
perform the forward procedure for Lyapunov exponents, 
besides of numerically solving the dynamical equations, 
one also needs to compute QR decompositions. For high- 
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Benettin steps, store $ and 
trajectory points 
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Iterate with (R-^)"\ 
compute CLVs 
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Figure 4. Computation of covariant Lyapunov vectors (CLVs). a) Method of LU factorization (see Sec. lIVBj) . and orthogonal 
complement method of Wolfe and Samelson (see Sec. IIV C)l . b) Iterative method of Ginelli et al. (see Sec. IIV D|) . 



dimensional systems good results are obtained with al- 
gorithms based on Householder transformations [1^ [l^ . 
The backward steps may in addition require an interpo- 
lation of the stored trajectory to find a solution of vari- 
ational equations with variable time steps. Finally, each 
column of A^(t„) is the null space of a corresponding 
rectangular submatrix of P. One of the most reliable 
methods of computation of the null space is based on the 
SVD . The null vector is identified as a right singu- 
lar vector corresponding to the vanishing singular value. 
Above we discussed that in principle in the case of degen- 
eracy of Lyapunov exponents one can obtain more than 
one null vector for one column A~(t„). But exactly iden- 
tical Lyapunov exponents are unlikely to occur in numer- 
ical computations, and, hence, multiple null vectors can 
(practically) never appear. It means that among right 
singular vectors we always have a preferable candidate 
with the smallest singular value. 

Implementations of QR decomposition and SVD in 
Fortran can, for example, be found in the well-known 
LAPACK library [13 ■ For a C++ implementations we 
refer to the ALGLIB NET hbrary [H]. Also this library 
provides implementations for many other platforms, such 
as Delphi and VBA. 



C. Orthogonal complement method of Wolfe and 
Samelson 

One of two first methods for the efficient computation 
of covariant Lyapunov vectors was suggested by Wolfe 
and Samelson |11J. Just as the LU method, their ap- 
proach utilizes the local property of the covariant vectors 
determined by Eq. (|43| . This equation can be written for 
the jthe vector as 

j 

i=l 
rn 

f,=T.V>t»tr (58) 

As above, the time dependence is not explicitly shown. 
Equating Eqs. ([58|) and ([57| and multiplying them by 
we can find 

i 

n=l 

Now we substitute this a^t,- in Eq. and multiply 

the resulting equation by ipj^ . Taking into account that 
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{Vkl]) = ^kj^ we obtain: 




where pik = (ip^if^) are elements of the matrix P (IT71) . 

In principle, this equation allows one to compute 
and to find the covariant vectors via Eq. (|57p . But this 
straightforward approach is not efficient. To compute 
the jth covariant vector, the coefficients a^^ are required, 
where fc = 1, 2, . . . , j. These coefficients depend on pik = 
{(pf(p^), where i = j,j + 1, . . . ,m. So, we need m — j + 1 
last vectors and j first vectors ip~ . The total number 
is always m + 1 . 

The key idea of Wolfe and Samelson to avoid this ob- 
stacle utilizes the orthogonality of P [Til , [soj . One can 
obtain the needed subspace spanned by the last (m—j+l) 
vectors by taking the orthogonal complement to the sub- 
space of the first (j — 1) vectors. In more detail, columns 
of P are orthogonal to each other, i.e., X^I^i PikPin = ^km 
where 6kn = 1 ii k ^ n and otherwise. This sum can 
be split at i = j as follows: 

m j-i 

PikPin = Skn - y^PifcPm- (61) 

i=j i=l 

The sum at the left hand side of this equation includes 
elements from the last rows of P, while the sum at the 
right hand side consists of the elements of the first rows. 
So, the sum in parentheses in Eq. (|60[l can be substituted 
as: 




Thus, to compute j unknown coefficients a"^, where 
n < we have to solve a set of j linear homogeneous 
equations 

X] ( ^PikPin I a,Tj =0 (j = 1, 2, . . . , m, k< j). 

n=l \i=l / 

(63) 

(We remind the reader that = for n > j.) Equa- 
tion ([63)) was suggested by Wolfe and Samelson to com- 
pute A^. It does not depend on the last rows of P, so 
that one needs j first backward vectors and j — 1 first 
forward vectors to compute j first covariant vectors. 

Later the method of Wolfe and Samelson was modified 
by Pazo et al. [l^ using the standard approach of com- 
putation of the forward and backward Lyapunov vectors, 
based on QR factorizations and on the backward itera- 
tions with the transposed propagator (these ideas were 
discussed in Sec. 



Changing the order of sums in Eq. (|63ll . we can write 
it in matrix form as 

P(l : J - 1, 1 :i)Tp(l : j - 1, 1 : j)A-(l : j, j) = 0. (64) 

Compare this equation with Eq. ([55)1 . We can see that 
solutions of Eq. ([55)1 constitute a subset of solutions of 
Eq. (I64p . But because we need only one solution at each 
j, and because our LU method finds such solution, we 
can conclude that the LU method works in the same way 
as the Wolfe and Samelson method, avoiding redundant 
matrix multiplication. 



D. Backward iterations, method of Ginelli et al. 

Almost simultaneously with Wolfe and Samelson, 
Ginelli et al. [l^l suggested a method based on asymp- 
totic properties of covariant vectors ([M)l . The under- 
lying idea of this method was described in the end of 
Sec. nil) but it cannot be directly implemented. Assume 
that we have backward Lyapunov vectors at ti. Theo- 
retically we can initialize Vj{ti) satisfying Eq. ([TU)) . and 
start the backward iterations using J-~^ . But in practice, 
due to numerical noise all these vectors shall belong to 
'S'm(ii) \ '5'„-i(ii)j because this set has the largest mea- 
sure. Hence, these iterations can provide only 7^. Due 
to the same reasons the forward iterations converge to 
7]^. The same is also true for the adjoint propagator. 

The key idea of Ginelli et al. is to perform the itera- 
tions in the space of projections onto backward Lyapunov 
vectors $~ (t). For a set of vectors initialized according to 
Eq. (|19p . the matrix of projections onto $ (t) is upper 
triangular and the iterations converge in the backward 
time. As follows from Eq. ([50)1 . the backward iterations 
with T(t\^t2)~^ in the space of projections onto $ (i) 
are equivalent to backward iterations with the upper tri- 
angular matrix R'^ (t 1 , ^2 ) ~ ^ ■ This mapping preserves the 
triangular structure of the matrix of projections, and we 
can perform as many backward iterations as we need al- 
ways staying within subspaces S~ {t\)\S~_-^(t\) . In other 
words, any upper triangular matrix iterated backward in 
time with R'^(ti,t2)~^ converges to A~(t). Note that 
since the subspaces S*^ (t) are spanned by the first j back- 
ward Lyapunov vectors, we are allowed to compute only 
J first covariant vectors without computing the rest of 
them. 

In a similar way we can compute the first j adjoint co- 
variant vectors, using the forward-time asymptotic ()24p . 
We start the procedure moving backward in time with 
the transposed propagator and computing forward Lya- 
punov vectors as described in Sec. [ill The triangular ma- 
trices R^(ti,t2) have to be stored. Then we turn round 
and start forward iterations R^(i„, t„+i)~^B(<„) = 
B(t„+i)[C^(i„,t„+i)]-i that converge to B+(t). 

A practical implementation of the method of Ginelli 
et al. might be the following; see the illustration in 
Fig. dKb). First, we perform the procedure for Lyapunov 
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exponents including forward iterations with ^(ii, ^2) and 
QR factorizations. This stage is preliminary and it is 
finished at iA when we decide that the orthogonal ma- 
trices Q{t) have converged to the matrices of backward 
Lyapunov vectors ^^(t). Starting from ^a, we continue 
the procedure, but now all the matrices #^(t„) and 
R'^(i„, see Eq. (|?T|l . are stored. This stage con- 
tinues until tB- The length of this stage depends on the 
number of points where we later want to compute the 
covariant vectors. After te we still proceed with the pro- 
cedure, but store only R'^(i„, i„+i). This stage must be 
sufficiently long to provide the convergence of the subse- 
quent backward procedure and it finishes at tc- At this 
point we initialize a set of arbitrary vectors, for which 
the property is fulfilled. In fact we just generate a 
random upper triangular matrix A. Using the stored ma- 
trices TV^ {t„,tn+i), we perform the backward iterations 
on the interval from tc to ^b- 

R^itn,tn+iy^A{tn+l) = A(t„)C (t„ , t„+i ) , (65) 

where the diagonal matrix C(i„, tn+i)^^ contains column 
norms of A. If ic ^ is sufficiently large, A(t„) con- 
verges to A^(i„). Now we pass the stage from to tA 
computing the covariant Lyapunov vectors via Eq. (j43p 
and using them as we need. Note, that this procedure 
allows one to compute not only the whole set of m co- 
variant vectors, but also as many of them as we want. 

As we already mentioned above, the columns of A~ (i„) 
can also be considered as covariant Lyapunov vectors, so 
that in some cases it is enough to consider these vectors 
without computation of r(i„). In this case the matrices 
$^(t„) do not have to be stored. 

The algorithm of backward iterations can suffer from 
ill-conditioned R"^, which manifests itself if one computes 
many (i.e., not just a few first) covariant Lyapunov vec- 
tors for a system with strong contraction. Typically, 
high-dimensional chaotic dissipative systems have sev- 
eral positive Lyapunov exponents of moderate magnitude 
while negative exponents can have large absolute val- 
ues. Because logarithms of diagonal elements of R"^ are 
proportional to local Lyapunov exponents, they can be 
sufficiently small. So, if a lot of covariant vectors corre- 
sponding to negative Lyapunov exponents are computed, 
the diagonal elements of R"^ can become small, and the 
whole matrix R"^, whose determinant is the product of 
its diagonal elements, can potentially be ill-conditioned. 
In turn this can influence the accuracy of computations. 

To avoid or at least minimize this problem one should 
first try to decrease the interval between QR orthogonal- 
izations. Another, also almost obvious recommendation 
is not to employ Eq. (|65p as it is, but compute iterations 
implicitly. Note that the implicit method is preferable re- 
gardless of the presence of ill-conditioned R"^. Namely, 
nonzero elements of the ith column of A(t„) can be com- 
puted as a solution of equation 

R-^(l:z,l:i)A„(l:z,i) = A„+i(l:i,i), (66) 



where R'^(l is a top left submatrix of R'^ and 

A„(l:i,i) top fragment of the ith column of A(i„). 
Computed in this way A„( : ,i) then has to be normal- 
ized. We see that the ith column of A(i„) is influenced 
only by the submatrix R (1 : i, 1 : i) that remains well- 
conditioned until i is sufficiently small. It means that 
even if R"^ has some small diagonal elements, errors that 
they can produce are not spread along the whole spec- 
trum, but infiuence only minor covariant vectors from its 
right part. 

When the trajectory passes close to tangencies of in- 
variant manifolds of an attractor. A(t„) becomes ill- 
conditioned, i.e., small values can appear on its diagonal. 
Because A(i„) is used to compute A(t„_i), small values 
can accumulate and vanish due to the numerical under- 
fiow. Then the zeros will be preserved in the course of 
iterations even if the trajectory goes far from the tan- 
gency points. This false indication of an exact tangency 
can be cured by adding a small amount of noise to the 
diagonal elements. 

E. Comparison of the methods 

Computation of covariant vectors requires saving of in- 
termediate matrices. We estimate the amount of the re- 
quired memory for the "worst" case when the whole set 
of m covariant vectors is computed. Let A"ab be the 
number of trajectory points where we are going to com- 
pute covariant vectors, i.e., the number of steps in the 
stage AB in Fig. [H It is reasonable to assume that this 
value depends on to, A'ab = A'ab('t^)j where to is the di- 
mension of the phase space. Denote the number of steps 
in the transient stage BC by ATbc- The convergence of 
columns of matrices to their asymptotic form during the 
transient stage is exponential with rates equal to differ- 
ences between corresponding Lyapunov exponents [ll| . 
For extensive chaotic systems these differences are pro- 
portional to 1/m; thus, the convergence time is propor- 
tional to TO. Altogether, the length of the transient stage 
can be estimated as A'bc = fcec^^j where k^c is an em- 
pirical constant, which depends on the particular system 
under consideration. 

For the LU method. Sec. IIVBI and for the method 
of orthogonal complement. Sec. IIV C[ the estimates are 
identical. Namely, we need Kab matrices each of the 
size TO^, and ATab + ATbc trajectory vectors of the size 
771, see Fig. HKa). Hence, the total amount of memory 
(in bytes) is Slu = ("7^(A'ab(to) -l-fcBc) +mKABim))b, 
where b is the number of bytes required to store one real 
number. For large to we have 

Bw ^m^KABim)b. (67) 

For example if the dimension is m = 100 and we want to 
compute A'ab = 1000 covariant vectors using double pre- 
cision numbers, i.e., b = 8, we need -Blu ~ 76 megabytes. 

For the method of backward iterations. Sec. IIVDI we 
need to save A'ab + A'bc triangular matrices R"^. each 
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of the size (m^ + m)/2, and T^'ab matrices $~ of the 
size m^, see Fig. ID^b). The total amount of memory 
can be estimated as B^i ~ (to^ (S-R'ab (^r^) + k^c'^^) + 
?Ti(i^AB('Ti)-l-fcBC"^))6/2. Keeping only the leading terms 
for large m we obtain: 



B 



BI 



m^(3XAB("^) + kBC'm)b/2. 



(68) 



For the same numerical values as in the example for LU 
method and at fcec = 1 we obtain, though higher, but 
close estimate: Bbi ~ 118 megabytes. Note however, 
that the amount of memory for the transient stage grows 
with m as fcBC™'^&/2 for the backward iterations method, 
while for two other methods it grows as A:bc'7i^6. Hence, 
the efficient application of the backward iterations re- 
quires closer attention to the minimization of the tran- 
sient stage length, otherwise, one can easily exhaust the 
available memory. 

In principle, all methods may suffer from a shortage 
of memory. One possible way to handle this problem 
is to save intermediate data to binary files. The disad- 
vantage of this approach is deceleration of computations 
due to the slowness of file operations. Alternatively, see 
Ref . [Ill , instead of keeping all necessary matrices moving 
forward in time, one can periodically (and sufficiently sel- 
dom to fit in the available memory) save snapshots of the 
procedure for Lyapunov exponents (i.e., the trajectory 
points of the basic system together with corresponding 
matrices $^). Then, moving backward, one periodically 
uses these snapshots to recompute forward steps and ob- 
tain missing data. Of course, this approach also slows 
down the computations, now due to the recomputations. 
To choose the preferable way one has to compare the av- 
erage time for writing to file and subsequent reading of 
one matrix with the time needed to recompute it. The 
result of comparison depends on the particular computer 
system. Note also that using the method of backward it- 
erations one can reduce the memory consumption if only 
the angles between covariant vectors are needed. As we 
already mentioned in Sec. lIVBj the triangle matrices A~ 
are suitable for finding the angles, and hence, in this case 
one does not need to save matrices 

Let us estimate the computation speed of the meth- 
ods presented (the straightforward intersection of the 
Oseledec subspaces is not taken into account). If all 
the methods have enough memory to avoid either using 
files or performing recomputing, the backward iterations 
are the fastest. Local methods of LU factorization and 
orthogonal complement loose the race on the backward 
stage B-A, see Fig. [H Each iteration is simultaneously 
a time step and also a computation of the covariant vec- 
tors. The time steps for local methods are performed 
via the procedure for Lyapunov exponents and also some 
time is required to compute the covariant vectors. 



EXAMPLES 



A. System with constant Jacobian matrix 



Consider a system with a constant Jacobian matrix 




(69) 



Since J is time-independent and has real eigenvalues, the 
Lyapunov exponents for this system simply coincide with 
the magnitude of its eigenvalues, Ai_2,3 = 1, —1, —3. The 
corresponding eigenvectors are simultaneously the covari- 
ant Lyapunov vectors, and the eigenvectors of (— J""") are 
the adjoint covariant vectors: 

T/2 \ 
T/2 1 -v/172 • (70) 

vX^ / 

D = ©'^r = diag[yT72, ^yT/3, V%^]- The propagator 
reads: 

/e^ e-^(l-e2^) \ 
J'ih,t2)^TLT-' = e-- , (71) 

\0 e-3^(e2^ - 1) e-^^J 

where t = t2 - h, and L = diag[c^i'^, e'^^'^, e^^'^]. For- 
ward and backward Lyapunov vectors can be computed 
as eigenvectors of far-future and far-past operators, re- 
spectively, directly from Eqs. and (finding the 
limits one has to keep constant norms of vectors): 





/T/2 ^1/2 0^ 

= I -v/172 V172 

ly 

(72) 

Note, that in accordance with Eq. (|43| . the first back- 
ward vector {1, 0, 0} and the last forward Lyapunov vec- 
tor {0, 0, 1}, coincide with the first and the last covariant 
vectors, i.e., with eigenvectors of J. One can also check 
that the logarithms of eigenvalues of the limit operators, 
i.e., the Lyapunov exponents, indeed coincides with the 
magnitude of the eigenvalues of J. The matrix P, as 
defined by Eq. (|T7l) . reads: 



0/2 -1/2 1/2 
/T72 1/2 -1/2 
./T/2 ^jt/2j 



(73) 



To compute covariant vectors via the LU method, we 
have to find the matrix A~. As follows from Eq. ([55]) . the 
first column of this matrix is always {1, 0, 0} while for the 
other elements we have a^^l^f^ — 032/2 = 0, a/j^/\/2 — 
^23/2 + 0;^3/2 = 0, and a{^/\f2 -|- 033/2 — 033/2 — 0. For 
the matrix B"*", needed to compute the adjoint covariant 
vectors, we construct equations according to Eq. ([5S)l us- 
ing P^: &l2/^/2 + 622/V2 = 0, 613/72 + 623/%/2 = 0, 
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— 5i3/2 + &23/2 + hj,-il\f2 = 0. Both of these equation 
sets have to be solved with the additional requirement of 
unit column norms: 

/I \ /I -v/172 1/2 \ 

A" = V273 7172 , B+ = 072 -1/2 . 

Vo ^172/ \0 71/2/ 

(74) 

One can check that Eqs. and are fulfilled, i.e., 
r = *-A- and = *+B+. 

The method of Wolfe and Samelson does essentially the 
same job. Computing A~ we have to multiply subma- 
trices of P by the transposed submatrices and construct 
equations; see Eq. ([64]). Similarly one can get B"'" and 
verify that the results coincide with Eq. ([74| . 

For the method of Ginelli et al. we find R'^(ti,t2) = 
[*"]'^:F(ti, ^2)*"; see Eq. ([31]). Since the iterations (|65l) 
converge in backward time, consider R'^(ti, ^2)"^: 

/e-" (e"-e--)/y2 (e"^ - e^)/\/2\ 
^^(tiMY^ = c^ e^^ -e^ . 

\ e^^ / 

(75) 

As follows from Eq. (|65l) . at t — ^ 00 the column norms 
of ~R^(tiM)~^ have to grow as e Indeed, it can 

be checked that the column norms of this matrix are 
asymptotically dominated by the terms e""^. e"^, and e^'^, 
respectively. If we normalize columns to the unit, the 
elements of this matrix converge to A^, see Eq. ([7H) . 
i.e.. we again obtain the covariant vectors. 

B. Generalized Henon map 

As second example we consider a generalized three- 
dimensional Henon map (40j 

<+i = a - \xlf - hxl 

x'^+' = (76) 

r^n+l _ n 
■^3 ~ '''2 ■ 

For a = 1.76 and & = 0.1 this system generates a hyper- 
chaotic attractor with Lyapunov exponents Ai = 0.225, 
A2 = 0.188, and A3 -2.716. 

Figure [5] shows the chaotic attractor, where the 
color of the points corresponds to det[P(l : 2, 1 : 2)] (see 
Sec. lIVB)) . Dark (red) colors indicate locations of the at- 
tractor where (almost) tangent CLVs occur and the sub- 
matrix P(l : 7, 1 : j) with j ~ 2 is (almost) singular. 

VI. CONCLUSION 

We presented an extensive description of modern 
achievements of Lyapunov analysis. The Lyapunov ex- 
ponents, the forward and backward Lyapunov vectors as 



well as covariant Lyapunov vectors were discussed in de- 
tail. 

The systematic approach allowed us to reveal a sym- 
metry in the structure of the tangent space and to in- 
troduce the concept of adjoint covariant vectors. There 
are tangent linear propagators that can be characterized 
by left and right singular vectors. When the propagators 
are considered on asymptotically growing time intervals 
these singular vectors converge to backward and forward 
Lyapunov vectors. One can also define adjoint propaga- 
tors that are associated with the same singular vectors, 
but have reciprocal singular values. The backward and 
forward Lyapunov vectors can be used as frameworks for 
two sets of Oseledec subspaces and for two adjoint Os- 
eledec subspaces that are orthogonal to the Oseledec sub- 
spaces. The main feature of these subspaces is the covari- 
ance with the tangent dynamics: the propagator maps 
each Oseledec subspace onto the corresponding Oseledec 
subspace associated with the image point of the trajec- 
tory, and the adjoint propagator does the same with the 
adjoint subspaces. Within these subspaces one can find 
vectors with the same property of covariance. There are 
covariant Lyapunov vectors whose exponential growth 
under the action of the propagators is characterized by 
Lyapunov exponents, and there are also adjoint covariant 
Lyapunov vectors that grow under the action of adjoint 
propagators with Lyapunov exponents of opposite signs. 

The adjoint covariant vectors are not independent 
characteristic vectors, because in principle one can al- 
ways compute them using the original covariant Lya- 
punov vectors. Important are the norm-independent an- 
gles between corresponding covariant and adjoint vectors. 
They provide a compact representation of the informa- 
tion provided by covariant vectors. In particular, homo- 
clinic tangencies between stable and unstable manifolds 
(characteristic for non-hyperbolic chaos) are indicated by 
orthogonality of corresponding original and adjoint vec- 
tors. 

An important result of our detailed analysis is an effi- 
cient method for computing covariant Lyapunov vectors. 
The basic idea of the method is an optimized LU decom- 
position of the matrix P consisting of scalar products of 
forward and backward Lyapunov vectors. Our approach 
is very close to the method by Wolfe and Samelson [Tl| . 
but its advantages are a more transparent explanation, 
and the explicit formulation of the matrix P which is 
interesting by itself. Moreover our approach is slightly 
more efficient because we avoid some redundant compu- 
tations. 

Using the matrix P, we present a method for detecting 
non-hyperbolicity of chaotic dynamics without explicit 
computation of the covariant vectors. In brief, the viola- 
tion indicator is the singularity of a j x j submatrix of P, 
where j is the number of positive Lyapunov exponents. 
The chaotic dynamics is non-hyperbolic if moving along 
a trajectory we encounter nearly singular submatrices. 

In presence of degenerate Lyapunov exponents all 
types of Lyapunov vector are not unique. We provide an 
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Figure 5. Attractor of the generalized Henon map Eq. 
(Color figure online) 



76)) . Dark (red) colors indicate closeness to homoclinic tangencies. 



analysis of this case. As for the forward and backward 
Lyapunov vectors, the standard algorithms can be used 
without modifications. Selection of an orthogonal initial 
matrix eliminates the ambiguity. Starting from different 
seed matrices, we can obtain different sets of vectors, but 
any one of them is appropriate. Moreover, in practical 
computations the degeneracy of the Lyapunov exponents 
manifests itself very weakly, especially for systems of high 
dimension. Typically, due to numerical errors all com- 
puted exponents are distinct, and one cannot identify 
degenerate exponents just by examining the computed 
spectrum. The same is true for the covariant vectors. 
Theoretically the degeneracy of the Lyapunov exponents 
can result in multiple sets of covariant vectors, but in 
practice the computations can be organized in a such 
way that one always obtains a unique appropriate solu- 
tion regardless of the degeneracy. 
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Appendix: Pseudocode for the LU method 

Input: nclv, number of computed covariant Lya- 
punov vectors; nstore, number of trajectory points 
where the covariant vectors are computed; m, dimension 
of the tangent space; dt, time interval between orthog- 
onalizations (normally, a multiple of time discretization 
step); nspend_att, nspend_fwd, nspend_bkw, steps to 
converge to the attractor, forward and backward vectors, 
respectively. 

Subroutines: solve_bas(), solving of the basic sys- 
tem; solve_lin_f wd() , 

solve_lin_trp(), action of forward and transposed 
propagators, respectively (see Sec. |ll|; null_vect(), 
computing a null vector (in the case of multiple 
solutions, an arbitrary null- vector can be taken); 
orthogO, QR-orthogonalization (matrix R is aban- 
doned); transpose 0, transpose of a matrix; randomO, 
generate random matrix or vector; A . B, multiplication of 
matrices A and B. 

Result: Gamma, array of nstore matrices m by nclv, 
whose columns are the covariant Lyapunov vectors. 

BEGIN clv_lu 

// *** ARRIVE AT THE ATTRACTOR *** 
CREATE u[l:m]=random(l,m) 
u=solve_bas (u,dt*nspend_att) 
// PRELIMINARY STAGE *** 

CREATE Q[l:m] [1 : nclv] =randoin(l ,m, 1 , nclv) 
Q=orthog(Q) 
FOR 1=1 TO nspend_fwd 
Q=solve_lin_f wd(Q ,u , dt) 
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Q=orthog(q) 
u=solve_bas (u , dt ) 
NEXT i 

// *** STAGE A-B *** 
CREATE PhiMns [1 :nstore] [l:m] [l:nclv] 
CREATE traj [1 : nstore+nspend_bkw] [l:m] 
FOR i=l TO nstore 

Q=solve_lin_f wd(Q ,u,dt) 

Q=orthog(Q) 

u=solve_bas (u , dt ) 

traj [i] =u 

PhiMns [i]=q 
NEXT i 

// *** STAGE B-C *** 
FOR i=l TO nspend_bkw 

u=solve_bas (u , dt ) 

traj [nstore+i]=u 
NEXT i 

// *** STAGE C-B *** 
// Now we use one column less 
RECREATE Q[l:m] [1 :nclv-l] =random(l ,m, 1 ,nclv-l) 
Q=orthog(Q) 

//We leave this cycle at 

// the (nstore+l)th trajectory point! 

FOR i=nspend_bkw TO 2 STEP -1 



u=traj [nstore+i] 
Q=solve_lin_trp(Q,u,dt) 
Q=orthog(Q) 
NEXT i 

// *** STAGE B-A *** 
CREATE P[l:nclv-1] [l:nclv] 
CREATE Gamma [1 : nstore] [l:m] [l:nclv] 
CREATE a[l:nclv] 

//We come into this cycle being at 
// the (nstore+l)th point 
// and take traj [i+1] , but not traj [i] . 
FOR i=nstore TO 1 STEP -1 
u=traj [i+1] 

Q=solve_lin_trp (Q , u , dt ) 
Q=orthog(q) 

P=transpose (Q) . PhiMns [i] 

Gamma [i] [l:m] [1] =PhiMns [i] [l:m] [1] 

FOR j=2 TO nclv 

a[l: j]=null_vect(P[l: j-1] [l:j]) 
Gamma[i] [l:m] [j] =PhiMns [i] [l:m] [1: j] .a[l: j] 
NEXT j 
NEXT i 
END 
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